En este trabajo se analiza la estructura multicapa de la red de transporte de la ciudad de Madrid a partir de dos conjuntos de datos reales proporcionados:
Los objetivos concretos son:
bus, metro,
bicimad, parking).Toda la parte empírica se basa exclusivamente en los ficheros:
transporte_madrid_consolidado.csvviajes_madrid.csvNo se utilizan ejemplos externos (como StarWars) para que el análisis esté completamente centrado en los datos de Madrid.
library(tidyverse)
library(janitor)
library(sf)
library(dbscan)
library(FNN)
library(igraph)
library(scales)
library(knitr)
Se asume que los ficheros CSV están en el mismo directorio que este documento RMarkdown. Si no es así, ajustar las rutas.
ruta_transporte <- "transporte_madrid_consolidado.csv"
ruta_viajes <- "viajes_madrid.csv"
transporte_raw <- readr::read_delim(
ruta_transporte,
delim = ";",
col_types = cols()
)
viajes_raw <- readr::read_delim(
ruta_viajes,
delim = ";",
col_types = cols()
)
glimpse(transporte_raw)
## Rows: 14,188
## Columns: 6
## $ stop_name <chr> " …
## $ stop_lat <dbl> 40.37835, 40.38341, 40.36813, 40.37773, 40.40753, 40…
## $ stop_lon <dbl> -3.75458, -3.62481, -3.62283, -3.67311, -3.59176, -3…
## $ transport_mode <chr> "bicimad …
## $ stop_id <chr> "382 -Parque Eugenia de Montijo, 2 …
## $ distrito_asignado <dbl> 2807911, 2807918, 2807918, 2807913, 2807919, 2807911…
glimpse(viajes_raw)
## Rows: 8,064
## Columns: 11
## $ fecha <chr> "2024-01", "2024-01", "2024-01", "2024-01", "2024-01…
## $ name <chr> "Centro "…
## $ edad <chr> "0-25", "0-25", "0-25", "0-25", "0-25", "0-25", "0-2…
## $ sexo <chr> "hombre", "hombre", "hombre", "hombre", "mujer", "mu…
## $ numero_viajes <chr> "0", "1", "2", "2+", "0", "1", "2", "2+", "0", "1", …
## $ personas <dbl> 90541756, 18809207, 73248991, 187551534, 75357987, 1…
## $ Codigo <dbl> 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, …
## $ Provincia <chr> "Madrid ", "Madrid ", …
## $ distrito <chr> "2807901 ", "2807901 ", "2807901 ", "2807901 ", "280…
## $ zona_pernoctacion <dbl> 2807901, 2807901, 2807901, 2807901, 2807901, 2807901…
## $ poblacion <dbl> 140644, 140644, 140644, 140644, 140644, 140644, 1406…
transporte_madrid_consolidadotransporte <- transporte_raw %>%
clean_names() %>%
mutate(
transport_mode = str_trim(transport_mode),
transport_mode = str_extract(transport_mode, "\\w+")
)
table(transporte$transport_mode)
##
## bicimad bus metro parking
## 626 12430 1050 82
summary(select(transporte, stop_lat, stop_lon, distrito_asignado))
## stop_lat stop_lon distrito_asignado
## Min. :40.28 Min. :-3.875 Min. :2807901
## 1st Qu.:40.40 1st Qu.:-3.711 1st Qu.:2807905
## Median :40.43 Median :-3.689 Median :2807909
## Mean :40.43 Mean :-3.683 Mean :2807910
## 3rd Qu.:40.46 3rd Qu.:-3.656 3rd Qu.:2807915
## Max. :40.56 Max. :-3.447 Max. :2807921
Variables clave:
stop_name: nombre de la parada/estación.stop_lat, stop_lon: coordenadas
geográficas.transport_mode: modo de transporte (bus,
metro, bicimad, parking).stop_id: identificador de parada.distrito_asignado: código de distrito (INE).Esto define los nodos base del sistema de transporte.
viajes_madridlibrary(readr)
library(dplyr)
library(janitor)
library(stringr)
# Ajusta la ruta a tu fichero si hace falta
viajes <- read_delim("viajes_madrid.csv",
delim = ";",
locale = locale(decimal_mark = ","),
show_col_types = FALSE) %>%
clean_names()
# Comprobamos que personas es numérica y razonable
summary(viajes$personas)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4515 37486 110989 140469 186820 737342
Variables clave:
fecha: mes de referencia.name: nombre del distrito.edad, sexo: características
sociodemográficas.numero_viajes: categoría de nº de viajes (0, 1, 2,
2+).personas: número de personas en esa celda.distrito: código de distrito.poblacion: población total del distrito.Construimos un indicador sencillo de viajes medios por
persona en cada distrito. Primero pasamos la categoría
numero_viajes a un valor numérico aproximado y luego
calculamos una media ponderada por el número de personas.
La variable personas se encontraba codificada con formato decimal europeo, por lo que fue necesario convertirla explícitamente a formato numérico. Asimismo, la variable numero_viajes es categórica (0, 1, 2, 2+), por lo que se transformó en una aproximación numérica asignando el valor 3 a la categoría “2+”. El indicador final de viajes medios por persona se calculó como una media ponderada por el número de personas en cada categoría, obteniéndose valores en torno a 2 viajes diarios por persona, coherentes con patrones conocidos de movilidad urbana.
library(dplyr)
library(knitr)
viajes_distrito <- viajes %>%
mutate(
# Pasamos las categorías de nº de viajes a un valor numérico aproximado
viajes_cat = dplyr::case_when(
numero_viajes == "0" ~ 0,
numero_viajes == "1" ~ 1,
numero_viajes == "2" ~ 2,
numero_viajes == "2+" ~ 3,
TRUE ~ NA_real_
)
) %>%
group_by(distrito) %>%
summarise(
# Total de personas (peso) en el distrito
personas_totales = sum(personas, na.rm = TRUE),
# Total de viajes ponderados
viajes_totales = sum(viajes_cat * personas, na.rm = TRUE),
# Viajes medios por persona
viajes_medios_por_persona = viajes_totales / personas_totales,
.groups = "drop"
)
kable(
head(viajes_distrito),
digits = 3,
caption = "Viajes medios diarios por persona y distrito"
)
| distrito | personas_totales | viajes_totales | viajes_medios_por_persona |
|---|---|---|---|
| 2807901 | 51888948 | 110673624 | 2.133 |
| 2807902 | 51542198 | 108276834 | 2.101 |
| 2807903 | 39251902 | 78164164 | 1.991 |
| 2807904 | 49039609 | 98112207 | 2.001 |
| 2807905 | 48349024 | 95741527 | 1.980 |
| 2807906 | 53679153 | 110757611 | 2.063 |
summary(viajes_distrito$viajes_medios_por_persona)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.956 2.004 2.038 2.041 2.063 2.135
El indicador de viajes medios diarios por persona presenta valores muy homogéneos entre distritos, situándose en un rango estrecho comprendido aproximadamente entre 1.98 y 2.13 viajes diarios por persona. Esta baja dispersión sugiere que, a escala distrital, los patrones de movilidad cotidiana en Madrid son relativamente uniformes, independientemente del volumen total de población o del número absoluto de desplazamientos. En particular, los distritos con mayor intensidad de viajes no necesariamente coinciden con los más poblados, lo que indica que la demanda de movilidad responde más a la estructura funcional del territorio y a la conectividad del sistema de transporte que al tamaño demográfico. Este resultado refuerza la idea de que la versatilidad multimodal y el papel estructural de ciertas zonas en la red pueden desempeñar un papel más relevante que la mera intensidad de viajes a la hora de caracterizar nodos clave del sistema de transporte urbano.
La idea es pasar de paradas individuales a zonas de intercambio (clusters de paradas cercanas), y luego construir una capa por modo de transporte.
# Creamos objeto sf con coordenadas
transporte_sf <- transporte %>%
st_as_sf(coords = c("stop_lon", "stop_lat"), crs = 4326, remove = FALSE) %>%
st_transform(25830) # sistema métrico (UTM zona 30)
coords_m <- transporte_sf %>%
st_coordinates() %>%
as.matrix()
# Clustering DBSCAN: paradas a menos de 150 m se agrupan
set.seed(123)
db <- dbscan(coords_m, eps = 150, minPts = 2)
transporte_sf$cluster_id <- db$cluster
# Tratamos el ruido (cluster 0) como clusters individuales
transporte_sf <- transporte_sf %>%
mutate(
cluster_id = if_else(
cluster_id == 0,
max(cluster_id) + row_number(),
cluster_id
)
)
# Zona de intercambio = centroide del cluster
zonas <- transporte_sf %>%
st_drop_geometry() %>%
group_by(cluster_id) %>%
summarise(
stop_lat = mean(stop_lat),
stop_lon = mean(stop_lon),
# distrito asignado más frecuente en el cluster
distrito_asignado = as.integer(names(which.max(table(distrito_asignado)))),
modos_disponibles = paste(sort(unique(transport_mode)), collapse = ","),
.groups = "drop"
)
nrow(zonas)
## [1] 1558
head(zonas)
Cada cluster_id representa una zona de
intercambio multimodal.
Representamos la distribución espacial aproximada (proyección geográfica simple):
zonas_sf <- zonas %>%
st_as_sf(coords = c("stop_lon", "stop_lat"), crs = 4326)
plot(zonas_sf["distrito_asignado"], main = "Zonas de intercambio por distrito")
modos <- sort(unique(transporte_sf$transport_mode))
zonas_modo <- transporte_sf %>%
st_drop_geometry() %>%
select(cluster_id, transport_mode) %>%
distinct() %>%
mutate(presente = 1) %>%
pivot_wider(
names_from = transport_mode,
values_from = presente,
values_fill = 0
)
zonas_attrs <- zonas %>%
left_join(zonas_modo, by = "cluster_id") %>%
mutate(across(all_of(modos), ~replace_na(., 0)))
kable(head(zonas_attrs), caption = "Atributos básicos de las zonas por modo")
| cluster_id | stop_lat | stop_lon | distrito_asignado | modos_disponibles | bicimad | bus | metro | parking |
|---|---|---|---|---|---|---|---|---|
| 1 | 40.38368 | -3.624114 | 2807918 | bicimad,bus,metro | 1 | 1 | 1 | 0 |
| 2 | 40.38096 | -3.728060 | 2807911 | bicimad,bus,metro | 1 | 1 | 1 | 0 |
| 3 | 40.38626 | -3.719795 | 2807912 | bicimad,bus,metro | 1 | 1 | 1 | 0 |
| 4 | 40.37886 | -3.732537 | 2807911 | bicimad,bus | 1 | 1 | 0 | 0 |
| 5 | 40.40855 | -3.672909 | 2807903 | bicimad,bus | 1 | 1 | 0 | 0 |
| 6 | 40.40380 | -3.706547 | 2807901 | bicimad,metro | 1 | 0 | 1 | 0 |
# Forzamos el mismo tipo en ambas columnas de distrito
zonas_attrs <- zonas_attrs %>%
mutate(distrito_asignado = as.integer(distrito_asignado))
viajes_distrito <- viajes_distrito %>%
mutate(distrito = as.integer(distrito))
# Join seguro
zonas_attrs <- zonas_attrs %>%
left_join(
viajes_distrito,
by = c("distrito_asignado" = "distrito")
)
summary(zonas_attrs$viajes_medios_por_persona)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.956 2.004 2.038 2.041 2.063 2.135
En muchas zonas puede que no haya dato (si el distrito no aparece o
hay problemas de asignación), pero para el análisis de red esto no es un
problema; simplemente tendremos algunos NA.
Para cada modo:
# Centroides en coordenadas métricas
zonas_centroides <- zonas_sf %>%
st_transform(25830) %>%
st_coordinates() %>%
as_tibble() %>%
rename(x = X, y = Y) %>%
mutate(cluster_id = zonas$cluster_id)
build_edges_modo <- function(modo, k = 3) {
zonas_modo_ids <- zonas_attrs %>%
filter(.data[[modo]] == 1) %>%
pull(cluster_id)
if (length(zonas_modo_ids) <= 1) {
return(tibble(modo = character(), from = integer(),
to = integer(), dist = numeric(), w = numeric()))
}
coords_modo <- zonas_centroides %>%
filter(cluster_id %in% zonas_modo_ids) %>%
arrange(cluster_id)
mat <- as.matrix(coords_modo[, c("x", "y")])
k_eff <- min(k, nrow(mat) - 1)
knn <- FNN::get.knn(mat, k = k_eff)
from <- rep(coords_modo$cluster_id, each = k_eff)
to <- coords_modo$cluster_id[c(knn$nn.index)]
dist <- as.numeric(knn$nn.dist)
tibble(
modo = modo,
from = from,
to = to,
dist = dist,
w = 1 / (1 + dist)
) %>%
distinct() %>%
filter(from != to)
}
edges_intra <- map_df(modos, build_edges_modo)
kable(head(edges_intra), caption = "Ejemplo de aristas intra-capa")
| modo | from | to | dist | w |
|---|---|---|---|---|
| bicimad | 1 | 288 | 649.2619 | 0.0015378 |
| bicimad | 1 | 300 | 321.5842 | 0.0031000 |
| bicimad | 1 | 12582 | 476.5779 | 0.0020939 |
| bicimad | 2 | 283 | 506.9860 | 0.0019686 |
| bicimad | 2 | 927 | 343.8273 | 0.0029000 |
| bicimad | 3 | 29 | 398.1358 | 0.0025054 |
g_list)nodos_ids <- zonas_attrs$cluster_id
build_graph_modo <- function(modo) {
e_modo <- edges_intra %>%
filter(modo == !!modo) %>%
select(from, to, w)
g <- graph_from_data_frame(
d = e_modo,
directed = FALSE,
vertices = zonas_attrs %>%
mutate(name = cluster_id)
)
missing_nodes <- setdiff(nodos_ids, as.integer(V(g)$name))
if (length(missing_nodes) > 0) {
g <- g + vertices(as.character(missing_nodes))
}
g
}
g_list <- map(modos, build_graph_modo)
names(g_list) <- modos
g_list
## $bicimad
## IGRAPH 6bce25f UN-- 1558 1401 --
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bce25f (vertex names):
## [1] 1-- 288 1-- 300 1--12582 2-- 283 2-- 927 3-- 29 3-- 735
## [8] 3-- 248 4-- 372 4-- 6182 4--15480 5-- 583 5-- 778 5-- 626
## [15] 6-- 296 6-- 1121 6-- 813 7-- 239 7-- 327 7-- 370 8-- 209
## [22] 8-- 152 8-- 145 11-- 1182 11-- 813 11-- 6031 12-- 6031 12-- 115
## [29] 12-- 4280 13-- 47 13-- 174 13-- 9493 14-- 1256 14-- 621 14-- 814
## + ... omitted several edges
##
## $bus
## IGRAPH 6bd1215 UN-- 1558 4036 --
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd1215 (vertex names):
## [1] 1--1289 1--1026 1--1303 2-- 652 2-- 283 2-- 29 3-- 819 3-- 390
## [9] 3-- 844 4-- 767 4--1288 4-- 638 5-- 158 5-- 422 5-- 969 7-- 20
## [17] 7-- 18 7-- 18 8-- 997 8-- 785 8-- 464 9-- 734 9-- 365 9-- 588
## [25] 10--1121 10--1011 10-- 239 11-- 31 11--1118 11-- 33 12-- 32 12-- 365
## [33] 12-- 532 13-- 296 13-- 938 13--1255 14-- 854 14-- 599 14--1182 15-- 43
## + ... omitted several edges
##
## $metro
## IGRAPH 6bd1ffc UN-- 1558 682 --
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd1ffc (vertex names):
## [1] 1--230 1--229 1-- 2 2--239 2--239 2--248 3--372 3--112 3--187
## [10] 6--187 6--248 6--917 7--287 7-- 68 7--145 11--753 11--220 11-- 49
## [19] 12--235 12-- 47 12--174 13--249 13-- 82 13--423 17--814 17-- 59 17-- 57
## [28] 18-- 57 18--516 18--564 24-- 30 24--446 24--223 27--222 27--437 27--118
## [37] 28--348 28--523 28--243 30--170 30--145 30--573
## + ... omitted several edges
##
## $parking
## IGRAPH 6bd2c18 UN-- 1558 161 --
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd2c18 (vertex names):
## [1] 7-- 29 7--14152 7-- 766 8--11008 8-- 239 8-- 245
## [7] 27-- 1052 27-- 814 27-- 814 28-- 446 28-- 657 28-- 164
## [13] 29-- 573 29-- 766 29-- 941 41-- 423 41--11008 41-- 562
## [19] 29-- 55 55-- 423 55-- 241 56-- 939 56-- 534 56-- 678
## [25] 67-- 613 67--14152 67-- 630 71-- 630 71-- 814 71-- 241
## + ... omitted several edges
muxVizEn esta sección generamos la visualización clave: un gráfico 3D en el que se ven las capas (modos) apiladas, con los nodos alineados entre capas.
Para ello usamos la función plot_multiplex3D de
muxViz. Se asume que muxViz está instalado; si
no lo está, el código no fallará, pero mostrará un mensaje.
if (requireNamespace("muxViz", quietly = TRUE)) {
library(muxViz)
layer_colors <- c(
bicimad = "#1a9641",
bus = "#d7191c",
metro = "#2c7bb6",
parking = "#fdae61"
)
node_multimodalidad <- zonas_attrs %>%
mutate(
modos_activos = rowSums(across(all_of(modos)))
) %>%
pull(modos_activos)
node_size_values <- rescale(node_multimodalidad, to = c(0.5, 2))
muxViz::plot_multiplex3D(
g.list = g_list,
layer.colors = layer_colors[names(g_list)],
layer.labels = "auto", # <- clave: escalar, no vector
node.size.values = node_size_values,
node.size.scale = 1.5,
edge.size.scale = 0.5,
layout = "fr",
show.aggregate = TRUE,
show.nodeLabels = FALSE
)
} else {
message("El paquete 'muxViz' no está instalado. Instálalo con devtools::install_github('manlius/muxViz') para ver la visualización 3D.")
}
Dado que la función plot_multiplex3D de
muxViz presenta problemas en la instalación utilizada (bug
interno con layer.labels), se ha implementado un gráfico
multicapa alternativo usando tidygraph y
ggraph, que representa las capas apiladas y las conexiones
intra- e inter-capa de forma equivalente.
# Gráfico multicapa alternativo con tidygraph + ggraph
library(tidygraph)
library(ggraph)
# Orden de las capas y un índice para apilarlas verticalmente
layer_order <- modos
layer_index <- setNames(seq_along(layer_order), layer_order)
# Distancia vertical entre capas (en coordenadas proyectadas)
offset <- (max(zonas_centroides$y) - min(zonas_centroides$y)) * 1.2
# Nodos multilayer: mismo cluster_id repetido por capa, apilado en vertical
nodes_multi <- purrr::map_dfr(layer_order, function(m) {
zonas_centroides %>%
transmute(
node = paste(cluster_id, m, sep = "_"), # id único por nodo-capa
cluster_id = cluster_id,
layer = m,
layer_idx = layer_index[m],
x = x,
y = y + (layer_index[m] - 1) * offset
)
})
# Aristas intra-capa: usamos edges_intra y las asociamos a cada capa
edges_intra_multi <- edges_intra %>%
transmute(
from = paste(from, modo, sep = "_"),
to = paste(to, modo, sep = "_"),
type = "intra",
layer = modo
)
# Aristas inter-capa: unen el MISMO cluster_id entre modos presentes (verticales)
inter_edges <- purrr::map_dfr(zonas_attrs$cluster_id, function(cid) {
# modos presentes en esa zona
fila <- zonas_attrs[zonas_attrs$cluster_id == cid, modos, drop = FALSE]
modos_presentes <- modos[as.logical(unlist(fila))]
if (length(modos_presentes) < 2) return(tibble())
comb <- t(combn(modos_presentes, 2))
tibble(
from = paste(cid, comb[, 1], sep = "_"),
to = paste(cid, comb[, 2], sep = "_"),
type = "inter",
layer = "inter"
)
})
edges_all <- bind_rows(edges_intra_multi, inter_edges)
# Grafo multilayer
g_multi <- tbl_graph(nodes = nodes_multi, edges = edges_all, directed = FALSE)
# Colores de capa (igual que antes)
layer_colors <- c(
bicimad = "#1a9641",
bus = "#d7191c",
metro = "#2c7bb6",
parking = "#fdae61"
)
ggraph(g_multi, layout = "manual", x = x, y = y) +
geom_edge_link(aes(alpha = type, linetype = type), show.legend = TRUE) +
geom_node_point(aes(color = layer), size = 1.2) +
scale_edge_alpha_manual(values = c(intra = 0.3, inter = 0.8)) +
scale_edge_linetype_manual(values = c(intra = "solid", inter = "dashed")) +
scale_color_manual(values = layer_colors) +
theme_void() +
labs(
title = "Red multicapa de transporte de Madrid",
subtitle = "Capas (modos) apiladas verticalmente y conectadas por zonas de intercambio",
color = "Capa",
linetype = "Tipo de enlace",
alpha = "Tipo de enlace"
)
Este gráfico muestra claramente:
bus, metro,
bicimad, parking) como planos separados.# ===========================================
# Visualización 3D multicapa (inter-capa con scatter3d lines)
# ===========================================
if (requireNamespace("plotly", quietly = TRUE)) {
library(plotly)
library(dplyr)
library(purrr)
library(scales)
library(tidyr)
# ---- 1. Capas y eje Z
layer_order <- modos
layer_index <- setNames(seq_along(layer_order) - 1, layer_order) # 0,1,2,3
# ---- 2. Nodos multilayer: solo zonas donde ese modo existe
nodes_multi_3d <- map_dfr(layer_order, function(m) {
zonas_attrs %>%
filter(.data[[m]] == 1) %>% # solo zonas con ese modo
select(cluster_id) %>%
inner_join(zonas_centroides, by = "cluster_id") %>%
transmute(
node = paste(cluster_id, m, sep = "_"),
cluster_id = cluster_id,
layer = m,
x_raw = x,
y_raw = y
)
})
# Normalizamos X, Y y asignamos Z
nodes_multi_3d <- nodes_multi_3d %>%
mutate(
x = rescale(x_raw, to = c(0, 1)),
y = rescale(y_raw, to = c(0, 1)),
z = layer_index[layer]
)
# ---- 3. Multimodalidad por zona
node_multimodalidad <- zonas_attrs %>%
mutate(modos_activos = rowSums(across(all_of(modos)))) %>%
select(cluster_id, modos_activos)
# Clusters con al menos 2 modos
clusters_interes <- node_multimodalidad %>%
filter(modos_activos >= 2) %>%
pull(cluster_id)
# ---- 4. Aristas inter-capa: todos los pares de modos presentes
edges_inter_3d <- map_dfr(clusters_interes, function(cid) {
fila <- zonas_attrs[zonas_attrs$cluster_id == cid, modos, drop = FALSE]
modos_presentes <- modos[as.logical(unlist(fila))]
if (length(modos_presentes) < 2) return(tibble())
comb <- t(combn(modos_presentes, 2)) # todos los pares posibles
tibble(
from = paste(cid, comb[, 1], sep = "_"),
to = paste(cid, comb[, 2], sep = "_"),
cluster_id = cid
)
})
cat("Número de aristas inter-capa generadas:", nrow(edges_inter_3d), "\n")
# ---- 5. Posiciones de nodos para las aristas
pos_nodes <- nodes_multi_3d %>%
select(node, x, y, z, layer)
edges_plot <- edges_inter_3d %>%
left_join(pos_nodes, by = c("from" = "node")) %>%
rename(
x = x,
y = y,
z = z
) %>%
left_join(pos_nodes, by = c("to" = "node"), suffix = c("", "_to")) %>%
rename(
xend = x_to,
yend = y_to,
zend = z_to
) %>%
filter(
!is.na(x), !is.na(y), !is.na(z),
!is.na(xend), !is.na(yend), !is.na(zend)
) %>%
select(x, y, z, xend, yend, zend, cluster_id)
cat("Aristas con coordenadas válidas:", nrow(edges_plot), "\n")
# ---- 6. Convertimos aristas a formato "long" para scatter3d lines
# Para cada arista: (x,y,z) -> (xend,yend,zend) -> NA como separador
if (nrow(edges_plot) > 0) {
edges_long <- edges_plot %>%
mutate(seg_id = row_number()) %>%
mutate(
x = map2(x, xend, ~c(.x, .y, NA_real_)),
y = map2(y, yend, ~c(.x, .y, NA_real_)),
z = map2(z, zend, ~c(.x, .y, NA_real_))
) %>%
select(seg_id, x, y, z) %>%
unnest(c(x, y, z))
} else {
edges_long <- tibble(x = numeric(), y = numeric(), z = numeric())
}
cat("Puntos totales en edges_long:", nrow(edges_long), "\n")
# ---- 7. Añadimos multimodalidad y tamaño de nodo
nodes_multi_3d <- nodes_multi_3d %>%
left_join(node_multimodalidad, by = "cluster_id") %>%
mutate(
modos_activos = replace_na(modos_activos, 1),
size = rescale(modos_activos, to = c(4, 11))
)
# ---- 8. Colores por capa
layer_colors <- c(
bicimad = "#1a9641",
bus = "#d7191c",
metro = "#2c7bb6",
parking = "#fdae61"
)
# ---- 9. Gráfico 3D
p3d <- plot_ly()
# Nodos
p3d <- p3d %>%
add_trace(
data = nodes_multi_3d,
type = "scatter3d",
mode = "markers",
x = ~x, y = ~y, z = ~z,
color = ~layer,
colors = layer_colors,
marker = list(size = ~size),
hoverinfo = "text",
text = ~paste(
"Zona:", cluster_id,
"<br>Modo:", layer,
"<br>Modos activos (total cluster):", modos_activos
)
)
# Aristas inter-capa como líneas rojas gruesas
if (nrow(edges_long) > 0) {
p3d <- p3d %>%
add_trace(
data = edges_long,
type = "scatter3d",
mode = "lines",
x = ~x, y = ~y, z = ~z,
line = list(color = "black", width = 4),
showlegend = FALSE,
hoverinfo = "none"
)
} else {
message("edges_long está vacío: no se han podido dibujar aristas inter-capa.")
}
p3d <- p3d %>%
layout(
title = "Red multicapa de transporte de Madrid (3D, enlaces entre capas)",
scene = list(
xaxis = list(title = "X (normalizada)"),
yaxis = list(title = "Y (normalizada)"),
zaxis = list(
title = "Capa (modo de transporte)",
tickvals = layer_index,
ticktext = layer_order
),
aspectmode = "cube"
)
)
p3d
} else {
message("Instala 'plotly' con install.packages('plotly') para ver la visualización 3D.")
}
## Número de aristas inter-capa generadas: 728
## Aristas con coordenadas válidas: 728
## Puntos totales en edges_long: 2184
Interpretación de la visualización 3D multicapa
La visualización tridimensional obtenida permite observar de forma explícita la estructura multicapa del sistema de transporte urbano de Madrid, donde cada plano horizontal representa un modo de transporte (bicimad, bus, metro y parking) y las conexiones verticales indican zonas de intercambio multimodal.
Separación clara de capas y estructura espacial
En primer lugar, se aprecia una separación nítida entre las cuatro capas, lo que confirma que cada modo de transporte presenta una distribución espacial propia:
Esta separación valida la decisión metodológica de modelar el sistema como una red multicapa, ya que una red agregada ocultaría estas diferencias estructurales.
Uniones verticales como indicador de multimodalidad
Las líneas verticales que conectan nodos entre capas representan zonas donde coexisten varios modos de transporte. Estas conexiones no aparecen de forma homogénea, sino que se concentran en un subconjunto reducido de zonas, lo que pone de manifiesto una fuerte heterogeneidad en la multimodalidad del sistema.
En términos de teoría de redes:
Estas zonas actúan como nodos altamente versátiles, ya que participan simultáneamente en varias capas.
Son puntos críticos para la integración funcional del sistema, permitiendo el transbordo eficiente entre modos.
Su posición estructural las convierte en puntos de control del flujo global de movilidad.
Visualmente, estas zonas destacan como “columnas” que atraviesan varias capas, reforzando la idea de que la conectividad global del sistema depende de un número limitado de nodos clave.
Versatilidad y roles multinodales
La figura ilustra de manera intuitiva el concepto de versatilidad de nodo:
Este patrón es consistente con la literatura sobre redes multicapas, donde se observa que la importancia de un nodo no puede evaluarse correctamente desde una sola capa. La visualización 3D confirma que los nodos versátiles no solo existen, sino que son estructuralmente visibles cuando se representan explícitamente las capas y sus interdependencias.
Implicaciones para vulnerabilidad y resiliencia
Desde el punto de vista de la vulnerabilidad del sistema, la visualización sugiere una conclusión relevante:
Este comportamiento es coherente con los modelos de cascadas de fallo en redes interdependientes descritos en la literatura (Buldyrev et al., 2010; Duan et al., 2019), donde la interdependencia aumenta tanto la eficiencia como la fragilidad del sistema.
Valor añadido de la representación multicapa
Finalmente, esta visualización demuestra el valor añadido del enfoque multicapa frente a representaciones tradicionales:
En conjunto, el gráfico 3D no solo actúa como herramienta exploratoria, sino como una evidencia visual clara de que el sistema de transporte de Madrid funciona como una red interconectada de capas, cuya eficiencia y resiliencia dependen de la correcta gestión de sus nodos multimodales clave.
Calculamos grado y fuerza (suma de pesos) para cada zona en cada modo.
centralidad_por_capa <- map2_df(
.x = g_list,
.y = names(g_list),
.f = function(g, modo) {
tibble(
cluster_id = as.integer(V(g)$name),
modo = modo,
degree = degree(g),
strength = strength(g, weights = E(g)$w)
)
}
)
kable(head(centralidad_por_capa), caption = "Centralidad (grado y fuerza) por capa")
| cluster_id | modo | degree | strength |
|---|---|---|---|
| 1 | bicimad | 5 | 0.0096077 |
| 2 | bicimad | 4 | 0.0090672 |
| 3 | bicimad | 4 | 0.0078289 |
| 4 | bicimad | 6 | 0.0134323 |
| 5 | bicimad | 7 | 0.0140216 |
| 6 | bicimad | 6 | 0.0139194 |
Distribución del grado por modo:
centralidad_por_capa %>%
ggplot(aes(x = degree, fill = modo)) +
geom_histogram(alpha = 0.6, bins = 30, position = "identity") +
facet_wrap(~ modo, scales = "free_y") +
theme_minimal() +
labs(
title = "Distribución del grado por modo de transporte",
x = "Grado",
y = "Frecuencia"
)
Definimos una medida de versatilidad basada en la entropía de Shannon de la distribución de grado entre modos.
mat_grado <- centralidad_por_capa %>%
select(cluster_id, modo, degree) %>%
pivot_wider(
names_from = modo,
values_from = degree,
values_fill = 0
) %>%
arrange(cluster_id)
grado_mat <- as.matrix(mat_grado[, modos])
rownames(grado_mat) <- mat_grado$cluster_id
calc_versatilidad <- function(vec) {
if (sum(vec) == 0) return(NA_real_)
p <- vec / sum(vec)
p <- p[p > 0]
H <- -sum(p * log(p))
H / log(length(vec))
}
versatilidad <- apply(grado_mat, 1, calc_versatilidad)
versatilidad_df <- tibble(
cluster_id = as.integer(names(versatilidad)),
versatilidad_grado = as.numeric(versatilidad)
) %>%
left_join(zonas_attrs, by = "cluster_id")
summary(versatilidad_df$versatilidad_grado)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.1477 0.0000 0.9983
top_versatiles <- versatilidad_df %>%
arrange(desc(versatilidad_grado)) %>%
slice_head(n = 15) %>%
select(cluster_id, versatilidad_grado, modos_disponibles, distrito_asignado)
top_especializados <- versatilidad_df %>%
arrange(versatilidad_grado) %>%
slice_head(n = 15) %>%
select(cluster_id, versatilidad_grado, modos_disponibles, distrito_asignado)
kable(top_versatiles, caption = "Top 15 zonas más versátiles")
| cluster_id | versatilidad_grado | modos_disponibles | distrito_asignado |
|---|---|---|---|
| 164 | 0.9983119 | bicimad,bus,metro,parking | 2807907 |
| 446 | 0.9976183 | bicimad,bus,metro,parking | 2807904 |
| 101 | 0.9968869 | bicimad,bus,metro,parking | 2807907 |
| 562 | 0.9963891 | bicimad,bus,metro,parking | 2807920 |
| 41 | 0.9955380 | bicimad,bus,metro,parking | 2807907 |
| 56 | 0.9949672 | bicimad,bus,metro,parking | 2807904 |
| 283 | 0.9949672 | bicimad,bus,metro,parking | 2807903 |
| 239 | 0.9946468 | bicimad,bus,metro,parking | 2807901 |
| 814 | 0.9934804 | bicimad,bus,metro,parking | 2807904 |
| 71 | 0.9927376 | bicimad,bus,metro,parking | 2807904 |
| 339 | 0.9927376 | bicimad,bus,metro,parking | 2807911 |
| 534 | 0.9927376 | bicimad,bus,metro,parking | 2807908 |
| 7 | 0.9926388 | bicimad,bus,metro,parking | 2807901 |
| 251 | 0.9892423 | bicimad,bus,metro,parking | 2807916 |
| 116 | 0.9863237 | bicimad,bus,metro,parking | 2807907 |
kable(top_especializados, caption = "Top 15 zonas más especializadas")
| cluster_id | versatilidad_grado | modos_disponibles | distrito_asignado |
|---|---|---|---|
| 9 | 0 | bus | 2807903 |
| 10 | 0 | bus | 2807909 |
| 15 | 0 | bus | 2807921 |
| 17 | 0 | metro | 2807921 |
| 19 | 0 | bus | 2807921 |
| 20 | 0 | bus | 2807921 |
| 21 | 0 | bus | 2807921 |
| 22 | 0 | bus | 2807908 |
| 26 | 0 | bus | 2807911 |
| 31 | 0 | bus | 2807921 |
| 32 | 0 | bus | 2807921 |
| 33 | 0 | bus | 2807921 |
| 34 | 0 | bus | 2807911 |
| 36 | 0 | bus | 2807911 |
| 38 | 0 | bus | 2807920 |
Interpretación cualitativa (a redactar en la memoria):
versatilidad_df %>%
ggplot(aes(x = viajes_medios_por_persona, y = versatilidad_grado)) +
geom_point(alpha = 0.4) +
geom_smooth(method = "loess", se = FALSE) +
theme_minimal() +
labs(
title = "Versatilidad de las zonas vs viajes medios por persona",
x = "Viajes medios por persona (distrito)",
y = "Versatilidad de grado"
)
Comprobamos cuántas observaciones completas hay y calculamos la correlación solo si tiene sentido estadísticamente (al menos 3 observaciones completas):
datos_cor <- versatilidad_df %>%
filter(!is.na(versatilidad_grado),
!is.na(viajes_medios_por_persona))
n_obs <- nrow(datos_cor)
n_obs
## [1] 1558
if (n_obs >= 3) {
cor_test <- cor.test(
datos_cor$versatilidad_grado,
datos_cor$viajes_medios_por_persona
)
cor_test
} else {
cat("No hay suficientes observaciones completas para estimar de forma fiable la correlación entre versatilidad y viajes.
")
}
##
## Pearson's product-moment correlation
##
## data: datos_cor$versatilidad_grado and datos_cor$viajes_medios_por_persona
## t = 0.13627, df = 1556, p-value = 0.8916
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.04621554 0.05310757
## sample estimates:
## cor
## 0.003454533
Esto evita el error de “not enough finite observations” y hace el análisis más robusto.
En lugar de implementar un modelo muy complejo de cascadas, realizamos un análisis sencillo pero informativo de robustez estructural:
# Matrices de adyacencia por capa
adj_list <- lapply(g_list, function(g) {
as.matrix(as_adj(g, attr = "w", sparse = FALSE))
})
# Aseguramos mismo orden de nodos
orden_nodos <- order(as.integer(V(g_list[[1]])$name))
adj_list <- lapply(adj_list, function(m) m[orden_nodos, orden_nodos])
adj_agregada <- Reduce("+", adj_list)
g_agregado <- graph_from_adjacency_matrix(adj_agregada, mode = "undirected", weighted = TRUE)
comp <- components(g_agregado)
tamano_gcc_inicial <- max(comp$csize)
tamano_gcc_inicial
## [1] 1558
Analizamos qué ocurre si vamos eliminando nodos con mayor grado en la
capa metro y miramos el tamaño relativo de la componente
gigante en la red agregada.
g_metro <- g_list[["metro"]]
deg_metro <- degree(g_metro)
orden_fallo <- names(sort(deg_metro, decreasing = TRUE)) # orden de fallo: más conectados primero
fracciones <- seq(0, 0.5, by = 0.05)
resultado_robustez <- map_df(fracciones, function(f) {
n_remove <- ceiling(f * length(orden_fallo))
nodos_fuera <- orden_fallo[seq_len(n_remove)]
g_agregado_reducido <- delete_vertices(g_agregado, nodos_fuera)
comp_red <- components(g_agregado_reducido)
gcc_red <- if (length(comp_red$csize) == 0) 0 else max(comp_red$csize)
tibble(
frac_eliminada = f,
tamano_gcc_rel = gcc_red / tamano_gcc_inicial
)
})
resultado_robustez %>%
ggplot(aes(x = frac_eliminada, y = tamano_gcc_rel)) +
geom_line() +
geom_point() +
theme_minimal() +
scale_y_continuous(labels = percent_format(accuracy = 1)) +
scale_x_continuous(labels = percent_format(accuracy = 1)) +
labs(
title = "Robustez de la red agregada ante fallos en la capa METRO",
x = "Fracción de nodos eliminados en METRO",
y = "Tamaño relativo de la componente gigante agregada"
)
Una caída rápida del tamaño de la componente gigante indica
alta vulnerabilidad ante fallos dirigidos en nodos
centrales de la capa metro.
En esta memoria se ha construido y analizado una red multicapas de transporte para Madrid utilizando exclusivamente los datos proporcionados:
bus,
metro, bicimad y parking.muxViz
que hace explícita la estructura en capas.metro.A partir de aquí, se podrían ampliar los análisis:
El código está preparado para ejecutarse directamente (asumiendo que los CSV y los paquetes necesarios están disponibles) y centrado únicamente en los datos reales de Madrid, tal y como se pedía en el enunciado.